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ABSTRACT 

The ratio of the vertical velocity dispersion to radial one (ct z /ctr) of self-gravitating 
bodies in various disc potentials is investigated through two different numerical meth- 
ods (statistical compilation of two-body encounters and TV-body simulations). The 
velocity dispersion generated by two-body relaxation is considered. The ratio is given 
as a function of a disc potential parameter, k/Q, where k and are the epicycle and 
circular frequencies (the parameters k/Q = 1 and 2 correspond to Kepler rotation 
and solid-body rotation). For 1 < k/Q 1.5, the velocity dispersion increases keep- 
ing some anisotropy (<7 z /<7r ~ 0.5-0.7) if the amplitude of radial excursion is larger 
than tidal radius, while <7 z /ctr <C 1 for smaller amplitude. On the other hand, for 
1.5 k/£1 < 2.0, we found isotropic state (<7 z /ctr ~ 1) in the intermediate velocity 
regime, while anisotropic state {<j z /(Jr < 1) still exists for higher and lower veloc- 
ity regimes. The range of the intermediate velocity regime expands with k/Q. In the 
limit of solid-body rotation, the regime covers all over the velocity space. Thus, the 
velocity dispersion generally has two different anisotropic states for each disc poten- 
tial (1 < «/f2 < 2) and one isotropic state for 1.5 k/H < 2 where the individual 
states correspond to different amplitude of velocity disper sion, while in the limit of 
solid-body rotation (k/H, — 2.0), entire velocity space is covered by the isotropic state. 

Key words: celestial mechanics, stellar dynamics - Galaxy: solar neighbourhood - 
methods: numerical - galaxies: kinematics and dynamics 



1 INTRODUCTION 

Gravitational interactions between bodies in a disc poten- 
tial tend to increase velocity dispersion of the bodies as well 
as diffuse the bodies radially. The increased random en- 
ergy is transferred from the potential energy through the 
radial diffusion. This process is called 'disc heating' for 
stars and molecular clouds in the Galactic gravitational field 
(e.g. Spitzer & Schwarzshild 1953; Binney & Tremaine 1987; 
Lacey 1991), and 'viscous stirring' for planetesimals in the 
protoplanetary disc or ring particles around planets (e.g. 
Stewart & Wetherill 1988; Ida 1990). Hereafter we use 'disc 
heating'. In general, disc heating results in anisotropic ve- 
locity dispersion. The radial and the vertical components of 
the velocity dispersion (o\r and <r z ) evolve with keeping a 
certain 'equilibrium' ratio, which is not generally unity. 

Ida & Makino (1992) showed through iV-body simu- 
lations that a z /an ~ 0.45 for self-gravitating planetesi- 
mals in the solar (Keplerian) potential. Numerical simula- 
tions of disc stars perturbed by massive melocular clouds 
in the solar neighbourhood showed ct z /ct_r ~ 0.6 (Villumsen 
1985; Kokubo & Ida 1992). Observations of stars in the so- 



lar neighbourhood show consistent anisotropy (e.g., Wielen 
1977; Chen, Asiain, Figueras & Torra 1997). 

In our galaxy, collective effects such as transient density 
waves would play an important role in velocity dispersion of 
stars in the solar neighbourhood (e.g., Barbanis & Woltjer 
1967; Binney & Lacey 1988; Jenkins & Binney 1990). Never- 
theless, it is important to clarify velocity dispersion created 
by (non-collective) two-body relaxation in a disc potential, 
because the two-body relaxation is one of the most basic 
processes in a self-gravitating disc system. Even if the col- 
lective effects dominate the disc heating in our galaxy, it is 
important to understand a competitive process, the heating 
by two-body relaxation. Furthermore, the two-body relax- 
ation may domminate in the central region of our galaxy 
or in other galaxies. In the present paper, we are concerned 
with dynamics regulated by the two-body relaxation in an 
'idealized' disc system with the potential characterized by 
a parameter k/Q, where k and Q. are horizontal epicycle 
frequency (see Eqs. (^)) and circular frequency of the disc 
potential, which indicates radial dependence of a disc po- 
tential. 

It is suggested that the equilibrium ratio (o~ z /o~r) 
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depends on k/Q. Lacey(1984) and Ida, Kokubo & 
Makino(1993) (hereafter IKM93) analytically investigated 
the equilibrium value of o z /or as a function of k/Q. They 
assumed that the evolution of the velocity dispersion is de- 
scribed by the sum of many independent two-body scat- 
terings in a disc potential with various initial conditions. 
They adopted the epicycle approximation (e.g. Binney & 
Tremaine 1987) and calculated orbital changes with the im- 
pulse approximation neglecting the external disc potential. 
IKM93 also suggested that the anisotropy is produced by 
deceleration of horizontal velocity at close approach due 
to shear motion (see section 2.2). For the Keplerian po- 
tential (k/Q = 1) and the galactic potential in the solar 
neighbourhood(fc/r2 ~ 1.4), IKM93 obtained the value of 
o z /or that are consistent with iV-body simulations (and ob- 
servations), while Lacey(1984) show significantly large val- 
ues (Fig. hi). IKM93 claimed that choice of the maximum 
impact parameter in the two-body formulae affects the equi- 
librium value of <t z /(Tr. IKM93 carefully chose the maximum 
impact parameter comparing with numerical orbital integra- 
tion to obtain smaller o z /or than that of Lacey(1984). How- 
ever, IKM93 obtained o z /or 7^ 1 in the limit of k/Q = 2 
(solid-body rotation). It would be reasonable to consider 
o z /or = 1 in the solid-body rotation, since shear motion 
vanishes. Lacey(1984) gave an argument with Jeans theo- 
rem to show o z j(TR = 1 in that case. Lacey(1984) obtained 
o'z/o'r = 1 for k/Q = 2, which is consistent with the above 
argument, although he failed to reproduce consistent values 
of Oz/or for smaller k/Q.. Things have been obscured be- 
cause of the lack of both observation and numerical work in 
the limiting case, k/Q ~ 2. 

To address this problem, we performed numerical sim- 
ulations in disc potentials with wide range of k/Q up to 
~ 2. In section 2, we calculate the disc heating as the sum 
of many independent two-body scatterings in a disc poten- 
tial, as Lacey(1984) and IKM93 did, but orbital changes 
are obtained by numerical orbital integrations. Our numer- 
ical calculations show an isotropic-dispersion (o z /or ~ 1) 
regime in velocity space if k/Q <; 1.5. We also found that 
the range of the regime expands with k/Q and the regime 
dominates all over the velocity space in the limiting case of 
the solid-body rotation (k/Q = 2). 

In section 3, to confirm our results obtained in section 2, 
we performed iV-body simulations of particles in various disc 
potentials. The results in section 2 and 3 are in good agree- 
ment with each other. In section 4, we summerize our results. 



2 VELOCITY EVOLUTION DUE TO MANY 
TWO-BODY SCATTERINGS IN DISC 
POTENTIALS 

2.1 Basic Formulation 

In this section, we consider a swarm of test bodies (parti- 
cle 1) gravitationally perturbed by field bodies (particle 2) 
in a disc potential field. We evaluate change rate of the ve- 
locity dispersion of the test bodies by statistically compiling 
the velocity changes in individual two-body encounters with 
field bodies, which are calculated numerically, following the 
method adopted by Ida (1990) and Kokubo & Ida (1992). 
Here we are also concerned with the cases with k/Q close to 
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Figure 1. The equilibrium ratio ctz/ctr as a function of k/Q. 
obtained by Lacey(1984) (solid line) and IKM93 (dotted line). 
Since the result of IKM93 is weakly dependent on the disc scale 
height owing to their choice of In A, we plot the case in which 
the disc scale height is larger than the tidal radius of a particle 
where their assumption is valid (the case where (i 2 ) 1 / 2 = 5, see 
section 2). 



2 while Ida (1990) and Kokubo & Ida (1992) only studied 
the cases with k/Q — 1.0 and 1.39. 

We assume that 'background bodies' generating the disc 
potential is continuously distributed and they do not con- 
tribute to two-body scattering. 

We adopt the epicycle approximation (see e.g. Petit & 
Henon 1986; Binney & Tremaine 1987): the velocity disper- 
sion is sufficiently smaller than rotational velocity around 
galactic centre. We use the following rotating coordinates: 



x= (R-a)/r s , 
y — (acj) — aQt)/r s 
z = z/r g , 



(1) 



where (R, <f), z) are cylindrical coordinates centred at the 
bottom of the disc potential, Q is the circular frequency 
at R = a. The normalization factor r g is defined by 



G(mi + 7712) 
Q 2 



(2) 



where mi and m.2 are the masses of the particle 1 and 2 and 
G is the gravitational constant. The radius r g corresponds 
to tidal radius within a numerical factor of O(l) except for 
k/Q ~ 2 (see Eq. @) 

In the epicycle approximation, the unperturbed orbits 
are given by (Petit & Henon 1986; Binney & Tremaine 1987) 
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and 



ay = ejsm{^t-Tj), 

a Q k 

Vi = ~2 b J + 2e J- cos (q* _ T > 

Zj =ijCos(—t-Wj), 



(4) 
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where time is scaled by and k and v are epicycle fre- 

quency and frequency of vertical oscillation which are de- 
fined by 



" = ' 0^ 



3fi 2 =2fi (r—\ 



dRJ R =a 



(5) 



where 4 > (-R, z) is an axisymmetric disc potential. For con- 
venience, we introduced a parameter a which indicates the 
strength of shear motion: 



K 



(6) 



The quantities e, i, b, r, u), and A are the constants of inte- 
gration. Equations (g) and (^) represent the particle motion 
as a combination of planar epicycle and vertical oscillation 
around the guiding centre rotating in a non-inclined circu- 
lar orbit. The quantities efi//t and r are the amplitude and 
the phase of the horizontal oscillation, respectively. Simi- 
larly, iQ/v and ui are those of the vertical oscillation. The 
first term of the right hand side of ijj represents the shear 
velocity, [Q(a + r g b) — fl(a)]a scaled by r g Sl (we assumed 
br g -C a). In the special case of Kepler rotation (k/Q = 1.0), 
the constants er g /a and ir g /a are called eccentricity and in- 
clination. 

From equations (Q), rms velocity averaged over an 
epicycle period of j-th particle is given by 



(JjR 



v h e 



(7) 



V2 



where ' - ' denotes time-averaging during the epicycle pe- 
riod. Then we define the velocity dispersion of a swarm of 
particle 1 (test bodies) as 



a R = -L= r 8 0, 
_ (i 2 ) 1 / 2 



_ r g Q, 
V2 



(8) 



,1/2 



where '( )' denotes ensemble-averaging (i.e. or = (<rf R ) 
and o z = <cr 2 z ) 1/2 ) We also call (e 2 ) 1/2 and (i?) 1/2 'nor- 
malized velocity dispersion' and pursue evolution of them. 
Kokubo & Ida (1992) considered gravitational scatterings of 
test bodies (disc stars) by many massive field bodies (giant 
molecular clouds) which are in non-inclined circular orbit. 
According to them, the evolution of the normalized velocity 
dispersion of the test bodies are given by 



d(ef) 2 
— — — n S 2rAl 
dt 

d<if> 



fi (ei , ii) Pheat (ei , ii) de 2 di\ , 



dt 



(9) 



n s2 r g Q / /i(ei,ii)<5heat(ei,ii)de 1 di 1 



where 

Pheat (e,i) 
Qhoat(e,i) 



-|b| 



A/ 



(2tt)2 
.2 a i, i drdaj 



db: 



db : 



p(e,i,b)db, 
q(e, i, b)db 



(10) 



2 1 1 (2tt) 2 

(we assumed b, r, and to are distributed randomly). The 



quantities p(e, i, b) and q(e, i, b) are introduced for later con- 
venience and written as 



p(e,i,b) 
q(e,i,b) 



, 2 a , T , drda; 

Ae 2' & W' 
A .2«,,i drdw 

Al 2 |6| (2^' 



(11) 



In Eqs. (^), n S 2 is the surface number density of the massive 
bodies and /i is the velocity distribution of the test bodies, 
which is normalized as 



/i(e, z)de 2 di 2 = 1. 



(12) 



Equations (|9j) are valid when 7712 2> mi . In more general case 
where the mass of a test body (particle 1) is comparable to 
that of a field body (particle 2), Eqs. ^ are revised as 



djef) 

dt 
d(ij) 

dt 



: n S 2r g Q ^ 
n 3 2r g Q. ^ 



m 2 



mi + 7T12 

m 2 
mi + (7i2 



(-Pheat) , 
(Qhcat) , 



(13) 



where 



(-Pheat) = / /(e,i)P hea t(e,i)de 2 di 2 



(Qhcat) = / /(e,i)<3heat(e,i)de di 



(14) 



(Ohtsuki 1998, Stewart & Ida 1998). In Eqs. ([jj), e and i 
are velocity dispersion of the relative motion defined by 



e cos t = e 2 cos r 2 — ei cos n , 
e sin r = e2 sin T2 — ei sin n , 
i cos w = 12 cos W2 — ii cos Wi , 
i sin w = 12 sin UJ2 — ii sin oji . 

Then, we obtain 



( e 2 ) = ( e 2 ) + ( e 2 ), 
(i 2 ) = (if) + (tl>. 



(15) 



(16) 



Under an assumption that both /1 and f 2 (velocity distri- 
butions of particles 1 and 2) are Rayleigh distribution, / is 
again Rayleigh distribution (Ohtsuki 1998, Stewart & Ida 
1998): 



/(e,z)de 2 di 2 



4ei 



exp 



■ exp 




Equations (|l^) and (Q imply that (Pheat) and (Qhcat), 
which are calculated only by orbital change in the relative 
motion, determine the ratio of the velocity dis per sion. It 
should be noted that integration with i in Eqs. (H) is also 
equivalent to that with v z and z, i.e., averaging with vertical 
velocity and height from disk mid-plane (Lissauer & Stewart 
1993). Strictly speaking, Eq. @ should include the terms 
expressing recoil of dynamical friction, which is proportional 
to (mi(e 2 ) - m 2 (ejj)) or (mi(i 2 ) - m 2 (il)) (Ohtsuki 1998; 
Stewart & Ida 1998), however, we neglect it, assuming the 
energy equipartition between particles 1 and 2 is already 
realized. In this case, it is shown that the recoil terms are 
much smaller than the heating terms (right hand sides of 
Eqs. (13)) if mi <C m 2 (Ida 1990). Furthermore, in the iden- 
tical particle case, Eq.(|13) is exact. 



© 0000 RAS, MNRAS 000, 000-000 



4 K. Shiidsuka and S. Ida 



The equations of the relative motion are given by (e.g. 
Icke 1982; Petit & Henon 1986; Kokubo & Ida 1992) 



2y 
■2x 



-x/r , 
-z/r 3 , 



(18) 



£1,2/2-2/1,22- 

2U/2 The last 



where time is scaled by 57 ~ , (x,y,z) — (x2 

zi), and r is scaled distance, r — (x 2 + y 2 -, ._ 

terms on the right hand side of Eqs. (hs|) represent the grav- 
itational interaction between two particles. Since we scale 
length and time by r g and 57, G(mi + m2)x/r' i is reduced 
to x/r 3 in these equations. The terms —2y and 2x are the 
Coriolis force. In the case of ft/57 = 1.0 and = 1.0, Eqs. 
( p^| ) are called Hill's equations which describe motion in the 
Kepler potential. 

We will numerically integrate Eqs. ([l8]) to evaluate 
the changes of orbital elements, in particular, e and i. 
The relative orbital elements are related to relative motion 
(x, y, z, x, y, z) in a similar way to Eqs. (0) and (^) as 



y = 



57 « 
o — e — cos( — t — r), 

ft 2 

a , „ 57 . , k , 
A - -bt + 2e— sm(-t-r), 



. 57 . . v 
i — sin — t 
v ' 



(19) 



and 



x = esin( — /: — r), 

a , „ 57 , ft , 

?/ = o + 2e — cos( — r — r), 

y 2 ^ ft v 57 ; ' 

i = i cos( — t — oj). 



(20) 



When the two particles are so far away that the mutual 
gravitational terms in Eqs. ( |l8| ) are negligible, the relative 
orbital elements are constants. When they approach each 
other, the orbital elements change through the mutual per- 
turbation. We show the examples of unperturbed motion, 
i.e., the motion with constant orbital elements in Fig. |^. In 
this figure, ft/57 = 1.87, e = 1.0, and b = 1.0 (solid line), 3.0 
(dotted line), and 5.0 (dashed line). 



2.2 Numerical Results of P^eat an d Qhcat 

We first investigate the behaviour of the 'elementary' quan- 
tities Phcat and Qhcat (Eqs. (|lo|)) since they show clearer 
physical properties than (Pheat) and (Qheat) and the aver- 
aging (Eqs. (|l^)) will not change the ratio of the velocity 
dispersion substantially. 

To obtain Pheat (Qheat), we calculate Ae 2 (Ai 2 ) with 
various sets of (b, r, uS) through orbital integration for each 
set of (e,i), following Ida (1990) and Kokubo & Ida (1992). 
As described in the previous subsection, only relative motion 
which obeys Eqs. ( |l8| ) is calculated. 

When relative distance r is large enough that mu- 
tual gravitation can be neglected, the orbital elements 
(e, i, b,T,u, A) are constants. We start our orbital integra- 
tion with sufficiently large y. A body approaches the other 
owing to shear motion. During the encounter, the orbital ele- 
ments are changed by mutual gravitational perturbation. We 
stop the integration when \y\ becomes large enough again. 
Changes of the orbital elements are determined as the differ- 
ence between orbital elements of initial and final states, i.e., 
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Figure 2. Examples of the unperturbed orbit on the x-y plane. 
For all cases, e = 1.0 and k/U = 1.87. Solid, dotted, and dashed 
lines are orbits with b = 1.0, b = 3.0, and b = 5.0, respectively. 



Ae 2 = e| llal — e 2 nitial . Since contribution in the integral ( |10[ ) 
from non-crossing distant encounters rapidly decreases with 
|6| in a disc system (Ida 1990; Hasegawa & Nakazawa 1990; 
Kokubo & Ida 1992) , initial b of orbits we calculated is con- 
fined in some finite regime. Furthermore, according to the 
symmetry of basic equations, the changes of e 2 and z 2 take 
the same values for the orbits with b and — b and those with 
u and uj + tv. Hence, we calculated orbits with < b < 6 max , 
< uj < n, and — n < r < n (for the value of 6 ma x, see 
below). 

The orbits are integrated with the fourth-order Her- 
mite scheme (Makino & Aarseth 1992). We also employed 
the algorithm developed by Emori, Ida, & Nakazawa (1993), 
where the part of deviation from the unperturbed epicycle 
orbit is numerically calculated while the part of the unper- 
turbed epicycle orbit is analytically calculated. 

We numerically calculated the heating rates with vari- 
ous K/fi: k/Q, = 1.00, 1.30, 1.58, 1.73, and 1.87. Since IKM93 
showed that anisotropy in the velocity dispersion is closely 
related to the shear motion, which depends on ft/51 but not 
on the parameter v/Q is fixed to 1.0. In the limit of 

ft/O — > 2.0, shear motion vanishes so that orbital integration 
becomes difficult. The case with k/Q ~ 1 is well investigated 
by several authors (e.g. Ida 1990, Kokubo & Ida 1992), hence 
we are mainly concerned with the results with k/Q. = 1.58, 
1.73, and 1.87. 

First we show the results for ft/57 = 1.87. In this case, 
we calculated Pheat and Qheat for 392 sets of e and i. For 
each set of (e, i), we integrated 10 4 - 10 6 orbits with different 
sets of (b, t, id). In Figs, ^-a to ^-c, we show obtained p(e, i, b) 
(solid lines) and q(e, i, b) (dotted lines) as functions of b. In- 
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tegration of p and q with b yields Phcat and Qhcat (Eqs. (Jic)) ) . 
In Figs. |-a, b, and c, (e,i) = (0.19,0.19), (1.23,1.23), and 
(10.1,5.04), respectively. Individual figures correspond to the 
results in the shear dominant region, the horseshoe domi- 
nant region, and the dispersion dominant region, which we 
define below. These figures suggest that only intermediate b 
contributes to the heating. The orbits with the intermediate 
b can closely approach each other. The orbits with smaller 
b turns back at distant y by the Coriolis force ('horseshoe 
orbits' (Brown 1911)), while those with larger b pass by at 
distant x. Ida (1990) and Kokubo & Ida (1992) showed that 
cumulative contribution (Phcat and Qheat) of distant encoun- 
ters with large b is negligible even if it is integrated over b 
to infinity. 

In the case of large e, the range of the strongly per- 
turbed orbits in the 6-space is extended by large amplitude 
of radial excursion, eQ/n, as shown in Fig. |j| In the unper- 
turbed orbits, the condition for orbit crossing is | b [< eQ/n 
(Eq.(|l^)). Hence we usually calculate orbits in the range of 
0.6 Si b 5s eVl/n + 2.5. In the cases of Figs. H-a, b, and c, 
we calculated 3 x 10 4 , 7 x 10 4 , and 4.3 x 1$ orbits. The 
numbers of calculated orbits are large enough that the in- 
tegrated values of Phcat (e,i) and Qheat(e,i) change by less 
than 5-10 per cent by the choice of calculated sets of initial 
conditions. 

In Fig. ^, we compile the calculated results for 392 sets 
of (e, i) with vectors. The angle between a vector and e-axis 
is determined by 

a , -l ( Qhcat/i 2 \ 

V Phcat /e 2 J 

while the length of the vectors are defined as a log[b(|P|/e 2 + 
\Q\/i 2 + 1)], where a and b are constants chosen for the 
vectors to be easily seen (in Fig. [l|, a = 0.2 and b = 200), 
and the factor 1 is introduced so that the length is always 
positive. These vectors show evolution trend of e and i. If 
a vector points to upper-right direction (8 — 45°), e and i 
increase keeping i/e constant. Since 

d(i/e) _ i , 1 di 2 1 de 2 
dt ~ 2e^P~dt ~ e 2 HP' 

i/e is kept constant when 

i 2 dt ~~ e 2 dt 
or equivalently, 

Qheat/i 2 = Phcat /e 2 . 

The evolution of e and i is divided into two steps. The first 
evolution is relatively rapid evolution toward the equilib- 
rium state of i/e, and the second one is the gradual increase 
in e 2 + i 2 keeping i/e constant. Fig. |^ suggests that the 
first step is faster than the second one. In the case where e 
and i are sufficiently large, we can quantitatively show the 
time-scale of the first step is much shorter than that of the 
second one (Appendix A). As a result, e and i would evolve 
along the line where vectors gather. Since de 2 /dt = CPheat 
and di 2 /dt — CQhcat, where C is some constant, we can 
integrate typical trajectory of a particle on the e-i plane 
using the data in Fig. ^ and their interpolation. The trajec- 
tory with initial condition (e,i) = (0.126,0.126) is plotted 
as dashed line in Fig. El Since the trajectory reaches the 



equilibrium state of i/e rapidly, it corresponds to a set of 
(e, i) in the 'equilibrium' state. 

Fig. ^ shows three different regions according to the 
manner of the evolution of e and i, in particular, equilib- 
rium ratio i/e: shear dominant region (e,i ,5s 0.4), horseshoe 
dominant region (0.4 S~ e,i Ss 2), and dispersion dominant 
region (e,i <; 2). In each region, the equilibrium ratio i/e 
is <C 1, ~ 1.1, and ~ 0.7, respectively, in the case where 
/t/O = 1.87. 

shear dominant region 

In the shear dominant region, the vertical and horizon- 
tal epicycle amplitude are small and shear motion domi- 
nates. Since shear motion is horizontal and orbits are bent 
before the orbits come close to each other, gravitational scat- 
tering takes place two-dimensionally. Accordingly, i hardly 
changes as shown in Fig. |§]-a. Dominant perturbation to 
e comes from orbits with b ~ 2. These orbits are distant 
'horseshoe'-type encounters as shown in Fig. ^-a. 

The orbital behaviour changes where the epicyclic os- 
cillation velocity, (e 2 + i 2 ) 1 / 2 ~ \Ple, is comparable to the 
shear velocity, namely, where ab/2 ~ y/2e (see Eqs. (pcj)). 
Hence the boundary should be e ~ 0.35, which is consistent 
with Fig. |. 

dispersion dominant region 

When e, i J; a, approach velocity is dominated by the 
random velocity v = (e 2 + i 2 ) 1//2 r g f2 rather than the shear 
velocity (abr g £l/2), and simultaneously, scattering occurs 
three dimensionally. Further, when the epicycle amplitude 
er g fl/n is larger than the tidal radius r t , where r t is defined 
by 

n - { cM? J (21) 

-1/3 

= r g a , 

that is, when e <; (k/S7)« -1//3 , orbits are not perturbed until 
the distance between bodies is much smaller than r t , since 
scattering cross section is small according to high relative 
velocity. Fig. |B|-c shows an example of orbital behaviour 
in this region. The orbit is hardly perturbed except when 
the distance is well smaller than r t . In this case, impulse 
approximation adopted by IKM93 is valid. They adopted 
Rutherford scattering formula neglecting a disc potential, 
and assumed incident motion to the two-body Rutherford 
scattering is described by the unperturbed motion given by 
Eqs. ^ and (|(i|). Actually, numerically obtained Pheat, 
Qheat, and e-i ratio in this region are in good agreement 
with those given by IKM93: in Fig. [|c, dashed line denotes 
p(e,i,b) and q(e,i,b) calculated by IKM93, and the inte- 
grated quantities, i.e., Phcat and Qhcat are in agreement with 
those obtained by IKM93 within accuracy of 10 per cent. In 
the region of small b, the analytical results deviate from the 
numerical ones, however, this deviation disappears when av- 
eraged over b. This is because some cancelation with regard 
to t or uj for fixed b in the analytical calculation would be 
transferred to cancelation for slightly different b by weak 
distant perturbation in the numerical calculation. 
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Figure 3. Dependence of p(e, i, b) (solid lines) and q(e, i, b) (dotted lines) on initial b. (a) Case with (e, i) = (0.19, 0.19). The b-dcpcndcnce 
is plotted every 0.05 (= 5b). To evaluate one point, 50 X 25 orbits with different r and u> are integrated, (b) Case with (e,i) = (1.23, 1.23). 
In this case, we varied 5b according to ((-dependence of p or q (we took smaller 5b where p or q rapidly changes): <56 = 0.01-0.1. (c) Case 
with (e, i) = (10.1, 5.04). As well as the case (b), we varied 5b according to b-dependence of p or q so that 5b = 0.01-0.2. For each b, we 
calculated 100 X 50 - 400 X 200 orbits with different r and u>. We also plot IKM93's result (dashed lines). 
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Figure 4. Directions and magnitude of the evolution of e and i on the e-i plane in the case where re/f2 = 1.87. The angles of the vectors 
are given by tan - 1 (Qhoat e2 /-Phcat* 2 ) an< i the length is max(a, a log[fe(|P|/e 2 + |Q|/i 2 + 1)]), where a = 0.2 and b = 200. Dashed line 
denotes the trajectory with initial condition (e,i) = (0.126,0.126). Solid line denotes the line i = 0.68e which is predicted by IKM93. 



In Fig.^, solid line denotes i/e = 0.68 which is predicted 
by IKM93. This is also consistent with numerical results for 
e £ (K/n)a- 1/3 ~ 2.4. 

In the dispersion dominant region, i/e in the equi- 
librium state is smaller than unity. The origin of this 
anisotropic velocity dispersion is explained as follows 
(IKM93). In the dispersion dominant region, relative mo- 
tion is approximated by unperturbed one as in Fig. |B|-c, and 
close encounter always takes place when the particle is mov- 
ing leftward (rightward) if b > (6 < 0). Hence at the 
moment of the closest approach, the horizontal component 
of the particle's incident velocity is always decelerated by 
the shear motion. On the other hand, the vertical motion is 
not affected by the shear and such a deceleration does not 
occur. At the moment of the closest approach, two-body 
scattering neglecting the tidal force can be applied, so that 
scattering tends to make (local) velocity isotropic, that is, 
equal energy is partitioned to each direction locally. Conse- 
quently, because of the deceleration of the horizontal motion 
at the closest approach, excessive energy is transferred to x- 
and y-directions compared to the z-direction. Hence in the 
dispersion dominant region, the equilibrium value of i/e is 
smaller than unity. 

As shown in Fig. 3-c, heating is dominated almost 
equally by orbits with 0.6 ^ b ^ efl/n. As described before, 
the upper limit eSl/n comes from the crossing condition of 
an unperturbed orbit. Hence, in large e case, orbits with 
correspondingly large b necessarily contribute to the heat- 
ing. In other words, encounters with large shear necessarily 



contribute. This is the case even in the limit of the solid- 
body rotation. Thus IKM93 obtained i/e ~ 0.8(< 1) even 
in the solid-body rotation case. As shown below, however, 
this anisotropic dispersion region practically disappears. 

horseshoe dominant region 

In the region with a/y/2 S~ e ^ (k/D.)^ 1 ^ 3 , approach 
velocity is dominated by random velocity as in the disper- 
sion dominant region. However, orbital behaviour is quite 
different from that in the dispersion dominant region. In this 
region, relative velocity is not as high as that in the disper- 
sion dominant region, so that relative motion is affected by 
distant perturbation similar to the shear dominant region. 
Figure 5-b shows an example of an orbit in this region. The 
orbit of guiding centre is gradually bent by Coriolis force 
until the orbit has b of different sign to turn back. Usually, 
such an orbit with e = i = is called 'horseshoe orbits'. Mo- 
tion of guiding centre is similar to 'horseshoe orbits' even in 
the case of e, i 7^ 0, since e and i are adiabatic invariants in 
distant region (Henon & Petit 1986). In this paper, we use 
the name 'horseshoe orbits' even if e, i ^ 0. 

As shown in Fig. 3-b, heating is dominated by orbits 
with 0.6 b ^ 2 in this velocity region. When b ^ 2, 
the horseshoe orbits are common, because Coriolis force 
dominates shear motion. However, in contrast with shear- 
dominated case, larger epicycle amplitude enables the bod- 
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Figure 5. Examples of orbit in the shear dominant region (a), the horseshoe dominant region (b), and the dispersion dominant region 
(c) in the case where n/Cl = 1.87. The solid lines denote the orbit, the dashed lines the guiding-center motion, and the circle with dotted 
line represents the tidal sphere (r = rt). Arrows represents the directions of the guiding-center motion. In each figure, initial orbital 
elements are (e,i,6) = (0.19,0.19,2.14) in (a), (e,i,6) = (1.23,1.23,2.14) in (b), and (e,i,b) = (10.1,5.04,4.2) in (c) (we omitted the 
other orbital elements such as phase variables). 
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10 r-T 
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e 

Figure 6. The same as Fig. 4 but k/Q = 1.58. The solid line 
denotes i = 0.60e predicted by IKM93 

ies to closely approach. In fact, close-encounting 'horseshoe 
orbits' as in Fig. 5-b contribute a lot to Ph ea t and Qheat- 

From the above orbital behaviour, the equilibrium ra- 
tio i/e ~ 1 is realized as follows. On the contrary to the 
dispersion dominant case, a close encounter always takes 
place when the particle is moving in the same direction as 
the guiding-centre motion (see Fig. ^-b). Hence, the parti- 
cle's horizontal motion is locally accelerated at the close en- 
counter. Consequently, the same argument as in the disper- 
sion dominant region predicts that i/e should became larger 
than unity. However, because the motion of the guiding cen- 
tre which is proportional to b is relatively slow, significant 
anisotropy does not appear. In the limit of the solid-body 
rotation, shear motion, and hence, the motion of the guiding 
centre slows down. Then i/e ~ 1 would be realized. 

So far we have considered the case where k/Q = 1.87. 
Our result in this case suggests that 

{the shear dominant region: e & a/v2, 

the horseshoe dominant region: . . 

a/V2&e£ [k/Q)^ 1 ' 3 , ( ' 
the dispersion dominant region: (k/S7)« _1//3 < e, 

where a — 4 — k 2 /Q 2 . To confirm the relation (0), we also 
calculated the cases with k/Q = 1.58 and 1.73 (Fig. | and 
Q). For k/Q = 1.58, we calculated 208 sets of (e,i), and for 
k/Q = 1.73, 218 sets of (e,i) are calculated. In each case, 
as in Fig. ^| we integrated typical trajectory of a particle on 
the e-i plane (dashed lines). The solid lines are also added 
as the result of IKM93. Expected boundaries of the regions 
in the individual cases are shown in Table |l|. Figs. ^, g and 
agree with Table |l[ 

The horseshoe dominant region shrinks as k/Q de- 
creases. The horseshoe dominant region would vanish for 
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Figure 7. The same as Fig. 4 but k/Q = 1.73. The solid line 
denotes i = 0.65e predicted by IKM93 



Table 1. Boundaries among three regions. 



k/Q 


a/V2 


(k/U)^ 1 / 3 


1.58 


1.0 


1.4 


1.73 


0.7 


1.7 


1.87 


0.35 


2.4 



a/V5 iS (ft/f^a -1 / 3 (equivalently, k/Q & 1.5). However, it 
should be noted that the 'horseshoe'-type close encounters 
still occur even if k/Q 1.5. We found when k/Q = 1.30, 
this effect makes the equilibrium ratio slightly larger than 
that predicted by IKM93 for e ~ 2-3. 

Since Ida (1990) and Kokubo & Ida (1992) only studied 
the cases with k/Q = 1.0 and k/Q = 1.39, respectively, they 
did not find the horseshoe dominant region. On the other 
hand, the horseshoe dominant region expands and domi- 
nates other two regions as k/Q approaches 2 (a — > 0). There- 
fore, in the limit k/Q — > 2, the region in which i/e ~ 1 is 
realized covers all over velocity space except for e — + oo. In 
the limit with e — > oo, IKM93's analysis would still be cor- 
rect. The IKM93's analysis, which well accounts for i/e in 
low k/Q case, is also valid in high k/Q case if e is sufficiently 
large, but the isotropic velocity region actually dominates in 
that case. Thus the contradiction stated in the introduction 
is solved. 

2.3 The Effect of Averaging on the Rayleigh 
Distribution Function 

Here we present the heating rates with the Rayleigh dis- 
tribution. This is necessary not only because the realistic 
velocity distribution must be considered, but also because 
we compare the results to those of iV-body simulations in 
the next section. 

Our numerical calculation of the heating rates Pheat and 
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Qheat is restricted within e, i & 12. By comparing numeri- 
cally obtained Pheat and Qhcat to those of IKM93, as stated 
in the last subsection, we find both are in agreement within 
accuracy of 10 per cent where e,i ii 7 for k/Q = 1.87. Ac- 
cordingly, in the region with ei>12oriJ>12, we use the 
analytical results of IKM93 in this case. In Fig. |^, we show 
the evolution diagram obtained from (Pheat) and (Qheat) 
in the case of k/Q = 1.87. The vectors drawn as a func- 
tion of the root mean squares, (e 2 ) 1//2 and (i 2 ) 1 / 2 in the 
same manner as in Fig. H (here, a — 0.2, b — 1000). Here, 
we use (e 2 ) 1 / 2 and (i 2 ) 1 ' 2 as normalized velocity disper- 
sion of particles and distinguish them from those of relative 
motion, (e 2 ) 1/2 and (i 2 ) 1/2 . Note that (e 2 ) 1/2 ((i 2 ) 1/2 ) and 
(e 2 ) 1/2 ((i 2 ) 1/2 ) are related by Eqs. (0). When mass of test 
particles 1 is much smaller than that of field particles 2, 
(e 2 ) 1/2 = (e 2 ) 1/2 and (i 2 ,) 1 ' 2 = (i 2 ) 1 ' 2 owing to energy 
equipartition. On the other hand, in the system of identi- 
cal mass, {el) 1 ' 2 = (e 2 )^ 2 /V2 and (i 2 ) 1 / 2 = (i 2 ) 1/2 /\/2. 
Hereafter we consider the system of identical mass in accor- 
dance with JV-body simulations in the next section. Hence it 
should be noted that e* is smaller than e by a factor V2- In 
this figure, we also plotted the typical trajectory in the same 
way as Fig. ^j. Although averaging smoothes the boundaries 
of the three regions observed in Fig. ^, the behaviour of the 
equilibrium ratio of the velocity dispersion is basically the 
same as Fig. |i| : In very small velocity case, (e 2 ) 1//2 3> (i 2 ) 1//2 , 
in intermediate velocity case ((e 2 ) 1 ^ 2 and (i 2 ) 1 ^ 2 are of or- 
der unity), the equilibrium value is almost unity, and in 
larger velocity case the equilibrium value is less than unity 
and seems to approach the value which predicted by IKM93 
(for example, in the equilibrium state, (i*) 1 ' 2 /(e») 1/ ' 2 ~ 1.0 
for (il) 1 ' 2 = 2.0, {e) 1/2 /{el) 1/2 ~ 0.95 for (i 2 ) 1 / 2 = 4.0, 
and (i 2 ) 1/2 /(e 2 ) 1/2 ~ 0.86 for (i 2 ) 1/2 = 9.0). It should be 
noted that the obtained equilibrium ratio seems to approach 
IKM93's much more moderately than the case of Fig. 0: In 
the case where the averaging is not done, our results (the 
equilibrium ratio) almost coincide with IKM93's for e J> 2, 
while with the averaging, our results do not coincide with 
IKM93's even when (e 2 ) 1/2 , (i 2 ) 1/2 ~ 10. In other words, 
the horseshoe dominant region influences even if (e 2 ) 1//2 and 
(i 2 ) 1 / 2 are much larger than unity as a result of the aver- 
aging on the velocity distribution. Actually, we found the 
influence of the horseshoe dominant region remains as long 
as (i 2 ) 1/2 £ 30 in the case of k/SI = 1.87 and (i 2 ) 1/2 <, 20 
for /t/fi = 1.58. On the other hand, in the case of Kepler 
rotation, where there is no horseshoe dominant region, we 
found that our numerical results coincide with IKM93's for 

(z 2 ) 1/2 ;>2. 



3 iV-BODY SIMULATIONS OF THE DISC 
HEATING 

3.1 Basic equations 

In section 2, we studied disc heating through statistical com- 
pilation of independent two-body scatterings. In this section, 
to check the results in section 2, we perform direct A-body 
simulations of particles in axisymmetric disc potentials. We 
consider N self-gravitating particles revolving in a disc po- 
tential. Then the motion of each particle is described as 
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Figure 8. Direction of the evolution of (e 2 ) 1 / 2 and (i 2 ) 1 / 2 plot- 
ted on the (e 2 ) 1 / 2 - (i 2 > 1/2 plane in the case of rc/fi = 1.87. The 
angle and length of each vector is determined in the same way as 
Fig. 4, but constants a = 0.2 and b = 1000. Similar to Figs. 4, 6, 
and 7, we plot the typical trajectory as the dashed line. 



At 



N 



(j = l,...,A0, 



(23) 



where the subscript j indicates particle's index, Vj , Xj , and 
rtij are the velocity vector, the position vector, and the mass 
of the particle j. The last term on the right hand side, Fj , is 
external force resulted from the disc potential. We consider 
the radial and vertical excursion of each particle are suffi- 
ciently small compared to its orbital radius in accordance 
with the study in the last section. In this case, the exter- 
nal force field Fj is expressed by two parameters of the disc 
potential, a and u/Q, as 



F vi = -afl 



F z i ~ —aQ 



(?) 
(?) 



a 2 \Qj a 2 



a 2 

2 ~2 



(24) 









\nj a 




a 2 \nj a 2 



where Cartesian coordinates (x, y, z) is centred at the bot- 
tom of the disc potential, R is the radial component of cylin- 
drical coordinates given by R — (x 2 + y 2 ) 1 ^ 2 . 

We integrate Eqs. (123) by using the fourth-order Her- 
mite scheme with the individual and hierarchical timestep 
(Makino 1991). The most expensive part of the Hermite 
scheme is the calculation of the force and its time-derivative 
whose cost increases in proportion to N 2 because we calcu- 
late the direct sum of all pairs. To reduce the computational 
time of this part, we employed a special purpose hardware, 
HARP-2 (Makino, Kokubo & Taiji 1993) and GRAPE-4 
(Makino, Taiji, Ebisuzaki & Sugimoto 1997). 
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Table 2. Initial parameters in TV-body simulations. 
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2000 4000 

Time [2%/ a] 

Figure 9. Time evolutions of (e 2 ) 1 / 2 (solid line) and (i 2 ) 1 / 2 
(dotted line) in the case of run 8. Time is scaled by the rotation 
period at R = a (2n/Q(a)). The dashed lines denote time evolu- 
tions of the normalized velocity dispersion predicted by d(e 2 )/di 
and d(i 2 )/di given by Eqs. (13) and the two-body results in sec- 
tion 2. 

3.2 Initial conditions of the swarm of particles 

We mainly investigated the k/Q dependence of the equi- 
librium ratio of velocity dispersion. The parameter v/Q is 
fixed to be unity in most cases, since IKM93 suggested that 
it does not affect the equilibrium ratio (we also did some 
iV-body simulations with v/Q 1 and found that the ra- 
tio is the same as in the case with v/Q = 1). We did 27 
iV-body simulations with nine different values of k/Q,. For 
each value of k/Q, we did several runs starting with different 
initial conditions. 

In Table ^, we summarize the initial conditions. We dis- 
tribute 1000 particles with identical mass randomly in the 
region a — Aa/2 < a < a + Aa/2. In most cases, the parti- 
cle masses are m — S x 10 _9 M g where M g is effective mass 
of the center defined by GM g /a 2 = Q 2 a. As suggested by 
the argument in section 2, the results would be indepen- 
dent of particle masses, if the results are scaled by r g . We 
took Aa < a to make simulation radially local, however, we 
took Aa sufficiently larger than rt (characteristic size of a 
particle's potential well) and epicycle amplitude for the edge 
effects to be negligible. We used enough number (N — 1000) 
of bodies that the bodies can closely approach each other. 
Initially, we set the same e, and i* (e,o, uo) for all particles, 
however, the velocity distribution converges to the Rayleigh 
distribution given by Eq. (^) in shorter time interval com- 
pared to the characteristic time-scale of the evolution of the 
velocity dispersion (i.e., two-body relaxation time T2b). 

3.3 Results of the iV-body simulations 

First we show the detailed results in two characteristic cases 
with k/Q — 1.30 and k/Q = 1.87. As stated before, orbital 
properties change near k/Q ~ 1.5. 

In Fig. |^, we show the time evolutions of (e 2 ) 1//2 (solid 
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k/CI 


u/a 
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Aa/a 
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1.00 


1.00 


0.6 


0.6 


0.29 
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1.00 


1.00 


0.6 
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0.27 
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0.27 
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0.13 
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0.27 
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0.6 


0.29 
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0.6 


0.29 
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0.27 


14 
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0.27 


15 
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1.00 


0.1 


0.1 


0.27 


16 
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1.00 


0,1 
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0.27 


17 
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1.00 


1.0 


0.5 


0.27 


18 


1.80 


1.00 
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1.0 


0.27 


19 


1.80 


1.00 


1.0 


2.0 


0.27 


20 


1.87 


1.00 


1.1 


0.6 


0.27 


21 


1.87 


1.00 


2.0 


1.0 


0.27 


22 


1.87 


1.00 


1.0 


2.0 


0.27 


23 


1.87 


2.00 


1.0 


2.0 


0.27 


24 


1.87 


4.00 


1.0 


2.0 


0.27 


25 


1.95 


1.00 


1.0 


0.5 


0.27 


26 


1.95 


1.00 


2.0 


1.0 


0.27 


27 


1.95 


1.00 


1.0 


2.0 


0.27 



line) and (i 2 } 1/2 (dotted line) of run 8 {k/Q = 1.30). The 
parameter k/Q = 1.30 corresponds to the galactic poten- 
tial in the solar neighbourhood. In the figure, time is scaled 
by the rotation period at a, namely, 2-7r/0(a). Although 
both {e 2 ) 1/2 and (i 2 } 1/2 keep growing, the ratio of (i 2 } 1/2 
to (e*) 1 / 2 seems to be constant in the later stage. The ratio 
(i 2 .) 1/2 /(e 2 ) 1/2 is plotted as a function of (e 2 +i 2 ) 1/2 in Fig. 
HQ, which is equivalently, time evolution of (i 2 ) 1//2 /(e 2 ) 1//2 . 
We also plot the results with different sets of (e*o,i*o) f° r 
the same k/Q (run 9 and 10). The ratio of the velocity dis- 
persions settles to a certain constant (~ 0.6) value inde- 
pendent of initial values of e* and i* after (e 2 + i 2 ) 1 ^ 2 ex- 
ceeds about 3. The equilibrium ratio gradually increases as 
(e 2 + i 2 ) 1//2 , which is consistent with the analytical argument 
in IKM93 that the equilibrium ratio has a weak dependence 
on (el+il) 1 ' 2 . 

The equilibrium ratio (i 2 ) 1 ^ 2 /(e 2 ) 1//2 ~ 0.6 is consis- 
tent with the observational value given by Wielen (1977) 
and Chen et al. (1996), statistical compilation of two-body 
encounters by Kokubo & Ida (1992), and iV-body simulation 
by Villumsen (1985). 

The time evolution of (e 2 ) 1/2 and (i 2 ) 1 '' 2 in the case 
of k/Q = 1.87 (run 22) is shown in Fig. [n]. Solid and 
dotted lines express (e 2 ) 1 ^ 2 and (i 2 ) 1 ^ 2 . The increase in 
(e 2 ) 1 / 2 is almost compensated by decrease in (i 2 ) 1 ^ 2 . The 
system evolve to the state of the equilibrium ratio, with- 
out increase in (e 2 + i 2 ) 1 ^ 2 , which is in contrast to the case 
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Figure 10. The evolution of ratio (i 2 } 1/2 /(e 2 } 1/2 against (e 2 + 
**} 1//2 in the cases of run 8, 9, and 10 (the same initial parameters 
but different sets of initial e*o and i*o)- This is equivalent to the 
time evolution, since (e 2 -M 2 ) 1 / 2 monotonically increases as time 
evolves. 
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Figure 12. The evolution of ratio (i 2 ) 1 / 2 /(e 2 } 1/ ' 2 against (e 2 
i 2 ) 1 / 2 in the cases of run 20, 21, and 22. 
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Figure 11. Time evolutions of (e 2 ) 1 / 2 (solid line) and (i 2 ) 1 / 2 
(dotted line) in the case of run 22. The dashed lines denote time 
evolutions of the normalized velocity dispersion predicted by the 
two-body results in section 2. 



of k/O = 1.30. Actually, the normalized random velocity 
(e 2 + i 2 ) 1 ' 2 only increases about 1 per cent throughout the 
simulation. This is consistent with the argument given in 
section 2 and appendix A (T ra tio T ran dom). In the poten- 
tial with k/Q, = 1.87, orbital motion is close to solid-body 
rotation. Since disc heating is caused by transformation from 
shear motion to random motion, the disc heating is weak in 
the potential with k/S1 = 1.87. Gravitational interactions 
mostly result in redistribution of random energy. 

In Fig. we show the evolution of (z 2 ) 1/2 /(e 2 ) 1/2 with 
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Figure 13. The time evolution of ratio (j 2 ) 1//2 /(e 2 ) 1 / 2 with 
various u/Q. Solid line corresponds to the results of run 22 
(u/n = 1.0). Dashed and dotted lines for run 23 (v/Sl = 2.0) 
and run 24 (u/Q = 4.0). 



different values of (e 2 + i 2 ) 1/2 (run 20, 21 and 22). This 
figure shows the equilibrium ratio settles to nearly unity 
independent of the initial conditions, as long as (e 2 ) 1//2 and 
(i 2 ) 1 ^ 2 are of order unity. 

We found that the distributions of the e* and i* evolve 
from the initial 5-function type function to the Rayleigh 
distribution in a time interval about < O.IT2B in both cases. 
Hence the assumption in Eq.(|l^) is valid. 

In Figs. [] and |ll], we also plotted the evolution of the 
normalized velocity dispersion calculated from d(e 2 )/dt and 
d(i 2 )/dt obtained in section 2 (Eqs. (|l3|)) as dashed lines. 
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Figure 14. Equilibrium value of (i 2 ) 1 / 2 / (e 2 ) 1 / 2 as a function of k/(1. Points with error bars are the results of the 27 TV-body simulations. 
Initial conditions of these simulations are presented in Table 2. Solid line corresponds to the isotropic state, i.e., (i 2 } 1//2 /(e 2 } 1//2 = 1.0, 
which is realized in the horseshoe dominant region. Dashed and dotted curves denote the results of IKM93 in the cases where (i 2 ) 1 / 2 = 1.5 
and (i 2 ) 1 / 2 = 5.0, respectively. Open squares are the equilibrium ratio obtained by the two-body results in section 2 at the mean velocity 
corresponding to the TV-body results. 



In both cases, the predicted results are in good agreement 
with the results of TV-body simulations. 

In addition to the runs with v/Q. — 1.0, we also did 4 
simulations with u/Q. = 2.0 and 4.0 (run 11, 12, 23, and 24). 
In the solar neighbourhood, v/£l ~ 4 (Binney & Tremaine 
1987). In these simulations, the other initial conditions are 
the same as run 8 for run 11 and 12, and the same as run 
22 for run 23 and 24. The time evolution of (i 2 ) 1/2 /(e*} 1/2 
in the case of k/Q, = 1.87 is shown in Fig. |l|. The solid line 
denotes the result of run 22 (u/Q. = 1.0), the dashed line 
is that for run 23 (v/Q. = 2.0), and the dotted line for run 
24 (v/Q. = 4.0). The variation of v/Q. does not affect the 
equilibrium ratio as expected, though the larger v/Q. results 
in the faster relaxation to the equilibrium state (This is be- 
cause the larger v/Q. leads to smaller disc scale height and 
therefore leads to more frequent scatterings among particles: 
see Eqs. (|^)). These results are the cases also for k/O, = 1.30. 

We also carried out TV-body simulations in the cases of 
k/Q = 1.0, 1.1, 1.2, 1.73, and 1.95, as summarized in Table 
§. In Fig. y, we plot the equilibrium value of {il) 1/2 / {e 2 ) 1/2 
obtained by the TV-body simulations with filled circles where 
error bars indicate standard deviation in time evolution. In 
the shear dominant region, we found timescale for heating 



(i.e. T2b) is comparable to or shorter than that for change 
of the ratio. Thus in the shear dominant region, there is 
no 'equilibrium' ratio. Accordingly, we are interested in the 
ratio in the dispersion dominant and the horseshoe dominant 
regions not in the shear dominant region. Since we cannot 
calculate the region with very large (e 2 ) 1 ^ 2 and (i 2 ) 1 ^ 2 for 
cpu limit, we only plot the results in the horseshoe dominant 
region in the cases of k/Q. <; 1.60. In the cases of k/Q. Ss 1.30, 
we plot the ratio in the dispersion dominant region, since 
there is no horseshoe dominant region. 

Together with the TV-body re- 

sults, the line (i 2 ) 1/2 /(e 2 ) 1/2 = 1 (the solid line) and the 
analytic lines in the dispersion dominant region by IKM93 
are plotted. The dashed and dotted lines are IKM93's results 
with (i 2 ) 1/2 = 1.5 and 5.0. In the plotted TV-body results, 
(i 2 ) 1/2 ~ 1.5-5 except run 15, 18, and 21. 

For k/Q. ~ 1.0-1.2, the TV-body results almost agree 
with the analytical line. For <; 1.87, the TV-body results 
show isotropy ((i 2 ) 1 / 2 / (e 2 ) 1 ^ 2 ~ 1), which is expected in the 
horseshoe dominant region. However, for 1.6 ^ re/Si ^ 1.8, 
the TV-body results are between the lines of IKM93 and the 
isotropy. As shown in section 2, the horseshoe dominant 
region is in limited velocity space in the case of relatively 
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Figure 15. Boundaries among three regions obtained from the relation (25) and corresponding equilibrium ratios. 



small k/Q. Since the velocity distribution (Rayleigh distri- 
bution) includes "dispersion-dominant" encounters more or 
less, its effect decreases the (i*} 1/ ' 2 /(e 2 } 1//2 from 1. On the 
other hand, the contamination of 'horseshoe'-type encounter 
would make (i 2 ) 1/2 /(e 2 ) 1/2 slightly larger than IKM93's re- 
sult as sug gested in section 2. To confirm this interpretation, 
in Fig. u4 we also plot the results (squares) of section 2 with 
Rayleigh distribution at the mean velocity corresponding to 
the iV-body results. They are consistent with the TV-body 
results. It is expected that numerical results would be con- 
sistent with IKM93's in sufficiently high mean velocity cases, 
because the effect of the horseshoe dominant region dimin- 
ishes then. In section 2, however, we suggested influence of 
the horseshoe dominant region remains when (i 2 ) 1/2 < 20- 
30 in the case where k/Q = 1.58-1.87, as a result of the 
averaging effect. Unfortunately, because of numerical diffi- 
culty we cannot directly examine such a high velocity re- 
gion with neither iV-body simulation nor the method in sec- 
tion 2. iV-body simulation with very high random velocity 
needs considerable cpu time since two-body relaxation time 
is proportional to the forth power of the random velocity 
in a disc system (Appendix A). The method in section 2 
is difficult to apply since scattering cross-section becomes 
very small so that we cannot obtain reliable results. In the 
limit of ft/fi — + 2.0, the effect of horseshoe dominant encoun- 
ters would always dominate and iV-body simulations always 
show (i 2 ) 1 / 2 /(e 2 ) 1 '' 2 ~ 1. 

In summary, the results of A^-body simulation are con- 
sistent with the results in section 2. Therefore, the physical 
argument in section 2 (mainly with quantities before veloc- 
ity averaging, P hea t and Qheat) would be valid. 



4 CONCLUSIONS 

We have studied the ratio of the velocity dispersion (<Tz/<tr) 
which is attained through mutual gravitational interaction 
among bodies in disc potentials. We examined the cases with 
Kepler rotation n/Q. = 1 to solid-body rotation n/Q. — 2, 



where k and are the epicycle and circular frequencies. 
Another parameter v/Q. does not affects the result, be- 
cause shear motion of particles, which is the origin of the 
anisotropic velocity dispersion, is independent of i//fi. We 
employed two different numerical methods, statistical com- 
pilation of two-body encounters (section 2) and iV-body sim- 
ulations (section 3). With the former method the physical 
properties are clearer and wider parameter range can be ex- 
amined, while the results are not direct and some assump- 
tions are introduced in the statistical compilation. On the 
other hand, the latter method is direct, although parameter 
range we can simulate is restricted by cpu time. The com- 
bination of the complementary methods would enable us to 
derive conclusive results. 

We found that the ratio becomes the equilibrium state 
much more quickly than the amplitude of the velocity dis- 
persion changes except when k/Q. is near 1. The equilibrium 
ratio depends on amplitude of velocity dispersion and disc 
potential parameter, k/£2. We found three characteristic ve- 
locity regimes: 



the shear dominant region: e <y a/V2, 

the horseshoe dominant region: 

a/V2<,e ^ (K/Q)a- 1/3 , 
the dispersion dominant region: (ft/f^a -1 / 3 JS e, 



(25) 



where a = 4— k 2 /Q, 2 and e corresponds to the amplitude of 
the random velocity of the relative motion of two particles 
normalized by r g Q. The velocity dispersion or and e are 



related as (e ) 



2,1/2 



2<Tfi/(r g f2) for a system of identical 



particle or (e 2 ) 1//2 = \/2 OR/(r s D.) for a system with large 
mass ratio such as stars and giant molecular clouds (see 
Eqs. 



and (|16[)). The characteristic radius r g is defined 
by Eq.(2|) and it is related to tidal radius r t as r g = a 1 ^ 3 ^. 

In the shear dominant region, shear motion dominates 
epicycle motion. Since shear motion is horizontal and or- 
bits are bent before the bodies come close to each other, 
gravitational scattering takes place two-dimensionally. As a 
result, only or is raised so that (t z /or <C 1. In the disper- 
sion dominant region, epicycle velocity is so large that orbits 
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are not perturbed until they closely approach each other. 
The close encounters are well approximated by Rutherford 
formula neglecting tidal force, as Lacey(1984) and IKM93 
did. As explained in section 2, energy equipartition in hor- 
izontal and vertical motion at the close encounters and the 
deceleration in horizontal velocity by shear motion lead to 
excessive or. IKM93 predicted <t z /(Tr = 0.5-0.8 {<j z /ur is 
smaller for smaller k/Q). Our numerical simulations agree 
with IKM93's prediction for relatively small fc/fi and sug- 
gest agreement even for larger k/Q. However, dispersion 
dominant region, where IKM93's prediction is valid, is over- 
whelmed by the newly found horseshoe dominant region 
in velocity space in the case of k/Q ~ 2 (See (^)). In 
the horseshoe dominant region, 'horseshoe'-type close en- 
counters dominates gravitational relaxation. In this case, 
o"z/o"j? ~ 1 is predicted. The physical reason for o z /<jr ~ 1 
is given in section 2. Due to the contamination from the 
encounters in other velocity regimes (the particles veloci- 
ties have Rayleigh distribution), iV-body simulations usually 
show a z /(jR smaller than 1. However, o z /<jr ~ 1 is actually 
shown in the case of k/CI ~ 2 as in Fig. [m] because the 
contamination diminishes as k/Q — > 2. 

In Fig. |l5| we draw a schematic figure on the ratio of 
velocity dispersion of self-gravitating particles in disc po- 
tentials from Kepler rotation to solid-body rotation. The 
boundaries of the three regions are given by the relation (p5j). 
In reality, the velocity distribution of particles makes the 
boundaries obscured through the averaging on the Rayleigh 
distribution. 
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APPENDIX A: 

We evaluate the time-scale for the equilibrium state of i/e 
to be realized and that for the random energy (e 2 + i 2 ) to 
be increased in the dispersion dominant region. The former 
is defined by 



Tratio — 



e d(i/e) 



At 

and the latter is 
1 



t random 



d(e 2 -M 2 



e 1 + i 2 



dt 



(Al) 



(A2) 



Accordingly, T ra tio and T ran dom are obtained from change 
rates of e 2 and i 2 . When e and i are large enough that the 
impulse approximation is valid, they are written as (Ida et 
al. 1993, Tanaka & Ida 1996): 



de^ 
dt 



= C 



di " -r 



3e 2 



K(X) 



3r 



+ V 



E(X) 



(A3) 



where K and E are the complete elliptic integral of the first 
and second kind, £ = 2fi/«, and A 2 = (1 - £~ 2 )e 2 /(e 2 + i 2 ). 
The factor C is defined by 



IT it 



ln(l + A 2 )- 



A 2 



- A 2 



1 



From Eqs. (|Al|), (|A2j), and (|A3|), We can rewrite T r 

^random 



-tratio — 
^random 



C 

e 2 + i 



C 



(2 + C)K(X)-3 



(1 + f 



2 i -2 

e + 1 



ATA)p, 



(A4) 
a tio and 

(A5) 



E{\) 



In above equations, when the non-dimensional factors mul- 
tiplied to (e 2 + i 2 )/C are order unity, it is easy to see both 
Tratio and Trandom are the same order as the Chandrasekhar's 
two body relaxation time. In the case where the disc po- 
tential is close to that of solid-body rotation (£ — > 1 or 
n/fi — * 2.0), however, Tandom tends toward infinity, while 
Tratio does not (we consider the case where particles have 
not reach the state of the equilibrium ratio yet). In Fig. 
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Figure Al. The ratio between T rat i and T random against k/Q. 
The solid line and the dashed line correspond to i/e = 1.0 and 
2.0. 



Al, r ra tio /T ran dom as a function of k/Q is plotted in the 
cases where i/e — 1.0 and 2.0 (note that T ra tio /Trandom is 
a function of i/e, but not e nor i). When k/Q is nearly 1 
(£ ~ 2), Tatio and Trandom are the same order. On the other 
hand, when k/Q is nearly 2 (f ~ 1), T ra ndom is much larger 
than Tatio- In the latter case, disc heating proceeds quasi- 
stationarily compared to the process to reach the state of 
the equilibrium ratio. 
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